Quantum Monte Carlo method using phase-free random walks 

with Slater determinants 
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We develop a quantum Monte Carlo method for many fermions that allows the use of any one- 
particle basis. It projects out the ground state by random walks in the space of Slater determinants. 
An approximate approach is formulated to control the phase problem with a trial wave function 
I^t). Using plane- wave basis and non-local pseudopotentials, we apply the method to Si atom, 
dimer, and 2, 16, 54 atom (216 electrons) bulk supercells. Single Slater determinant wave functions 
from density functional theory calculations were used as Y&t) with no additional optimization. The 
calculated binding energy of Si2 and cohesive energy of bulk Si are in excellent agreement with 
experiments and are comparable to the best existing theoretical results. 
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Quantum Monte Carlo (QMC) methods based on aux- 
iliary fields (AF) are used in areas spanning condensed 
matter physics, nuclear physics, and quantum chemistry. 
These methods jl], |[ allow essentially exact calculations 
of ground-state and finite-temperature equilibrium prop- 
erties of interacting many fermion systems. The required 
CPU time scales in principle as a power law with system 
size, and the methods have been applied to study a va- 
riety of problems including the Hubbard model, nuclear 
shell models, and molecular electronic structure. The 
central idea of these methods is to write the imaginary- 
time propagator of a many-body system with two-body 
interactions in terms of propagators for independent par- 
ticles interacting with external auxiliary fields. The in- 
dependent particle problems are solved for configurations 
of the AF and averaging over different AF configurations 
is then performed by Monte Carlo (MC) techniques. 

QMC methods with auxiliary fields have several ap- 
pealing features. For example, they allow one to choose 
any one-particle basis suitable for the problem, and to 
fully take advantage of well-established techniques to 
treat independent particles. Given the remarkable de- 
velopment and success of the latter |J, it is clearly very 
desirable to have a QMC method that can use exactly 
the same machinery and systematically include correla- 
tion effects by simply building stochastic ensembles of the 
independent particle solutions. Vigorous attempts have 
been made from several fields to explore this possibility 
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A significant hurdle exists, however: except for special 
cases (e.g., Hubbard), the two-body interactions will re- 
quire auxiliary fields that are complex. As a result, the 
single-particle orbitals become complex, and the MC av- 
eraging over AF configurations becomes an integration 
over complex variables in many dimensions. A phase 
problem thus occurs which ultimately defeats the alge- 
braic scaling of MC and makes the method scale expo- 
nentially. This is analogous to but more severe than the 
fermion sign problem with real AF f% Bfl or in real-space 



methods |l(J. No satisfactory, general approach exists 
to control the phase problem. As a result, only small 
systems or special forms of interactions can be treated. 

In this paper we address this problem. We develop 
a method for many-fermions that allows the use of any 
one-particle basis. It projects out the ground state by 
random walks in the space of Slater determinants. The 
phase problem is eliminated with an approximation that 
relies on a trial wave function \^t)- We demonstrate 
the method by applying it to electronic systems using a 
plane-wave basis and non-local pseudopotentials, which 
can be implemented straightforwardly in this method. 
We calculate the binding energy of Si2 and the cohesive 
energy of bulk Si using fee supercells consisting of up to 
54 atoms (216 electrons). These calculations represent 
the first application of AF-based QMC to solids. The re- 
sults are in excellent agreement with experiments and are 
comparable to the best existing theoretical results. Par- 
ticularly worth noting is that our results were obtained 
with a trial wave function which is a single Slater deter- 
minant formed by orbitals from density functional theory 
(DFT) calculations (with the local density approximation 
(LDA)), with no additional parameters or optimization. 

The Hamiltonian for any many-fermion system with 
two-body interactions can be written in any one-particle 
basis in the general form 
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H = H 1 +H 2 
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Vijkic\c]ckCi, (1) 



where N is the size of the chosen one-particle basis, and 
cj and Ci are the corresponding creation and annihilation 
operators. Both the one-body (Ty) and two-body matrix 
elements {Vijki) are known. 

To obtain the ground state \^g) of H, QMC methods 
use the imaginary time propagator e~ rH acting on a trial 
wave function |* T ): lim, woo (e~ ri ^)™ |* T ) oc |* G ). |* T ) 
must not be orthogonal to \^a), and we will assume that 
it is of the form of a single Slater determinant or a linear 
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combination of Slater determinants. The time step r is 
chosen to be small enough so that Hi and H2 in the 
propagator can be accurately separated with the Trotter 
decomposition. 

The propagator e~ rHl is the exponential of a one-body 
operator. A propagator of this form acting on a Slater 
determinant is straightforward to calculate, and it sim- 
ply yields another determinant. The two-body propaga- 
tor e~ rH2 can be expressed as an integral of propaga- 
tors of this form, as follows. Any two-body operator can 
be written as a quadratic form of one-body operators: 
H 2 = ^|Ea^«''ai where A Q is a real number and v a 
is a one-body operator. The Hubbard-Stratonovich (HS) 
transformation [TL ill then allows us to write 
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Introducing vector representations a = {<7i, 02, • • •} and 
v = {\/Ai v\, V^2i>2, ■ ■ •}, we have the desired form 



-tH 



J P{a) B(a) da, 
where -P(er) is the normal distribution in Eq. (^|) and 
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B(a) = e~ TH ^ e 



tHx/2 p ^a-w g-rffi/2 



(4) 



is a one-body propagator. 

The imaginary-time propagation thus requires evalu- 
ating the multidimensional integral in Eq. (^|) over time 
slices n and the corresponding auxiliary fields. MC tech- 
niques are the only way to evaluate such integrals effi- 
ciently. We use a random walk approach J9|. In each 
step, a walker |</>), which is a single Slater determinant, 
is propagated to a new position |</>'): \4>'{a)) = B{a)\4>), 
where a is a random variable sampled from P{a). After 
a sufficient number of steps (iterations), the ensemble of 
random walkers is a MC representation of the ground- 
state wave function: = ^ \cj)'). 

In general A Q cannot be made all positive in Eq. (j|) 
jn|. The one-body operators v are therefore complex. 
As the projection proceeds, the orbitals in the random 
walkers will become complex. As a result, the statistical 
fluctuations in the MC representation of \^g) increase 
exponentially with projection time (3 = nr. This is the 
phase problem referred to earlier. It is of the same origin 
as the sign problem that occurs when B(a) is real. The 
phase problem is more severe, however, because for each 
instead of a +\<fi) and —\(f>) symmetry ||, there is 
now an infinite set {e l6 \<j)}} (9 £ [0, 2n)) from which the 
random walk cannot distinguish. At large /3, the phase 
of each becomes random, and the MC representation 
of |*g) becomes dominated by noise. This problem is 
generic, and the same analysis would apply if we had 
chosen, instead of the random walk, the standard AF 



FIG. 1: Illustration of the phase problem and constraints to 
control it. The total valence energy (in Ry) of an fee Si prim- 
itive cell (2 atoms) is shown as a function of projection time 
/3 = nr, with r = 0.05 Ry -1 . Unless otherwise indicated, 
10, 000 walkers are used. Increasing the number of walkers 
from 50 to 10, 000 only slightly delays the onset of the phase 
problem. Simple generalization of the constraint that worked 
well for real determinants leads to poor results. The new 
method gives accurate results (note the agreement with the 
solid free projection curve, which is exact, until the latter 
becomes too noisy at /3 ~ 1.5). 



QMC sampling approach E| . In Fig. 0, the curves labeled 
"free projection" illustrate the phase problem. 

Existing fixed-node type approximations have often 
worked very well to control the sign/phase problem in 
real space |l3|, |l4| or in Slater determinant space when 
the propagator is real [§. The phase problem here is 
unique because not only do the determinants acquire 
overall phases, but the internal structures of their or- 
bitals become complex. The real-space analogy would 
be to have walkers whose coordinates become complex. 
This makes straightforward generalization of existing ap- 
proaches ineffective. For example, similar to the con- 
strained path approximation |^| we could impose the con- 
dition Ke(^T\4>) > 0. Or, in the spirit of the fixed-phase 
approximation in real space |l4[ ] we could project the 
walker by including a factor cos(A#) in the weight, where 
A8 is the phase of (^t\4>') I '(^t|0) ■ They give similar re- 
sults and do not work well ||. The former is shown in 
Fig. |l] ( "simple constraint" ) . Importance sampling with 
Re^Tl^) or K^tI^)! does not change the results Jl5[ . 

To formulate a new method that can better separate 
the overall phase from the determinant, we first borrow 
from the idea of importance sampling |l6[| , although our 
choice of the so-called importance function, (^t\4>): is 
actually complex. We modify Eq. (||) to obtain the fol- 
lowing new propagator for |0): 



(tf T |<£'(o- - 5-)) P{a ~ a)B(a - a)- 



da, (5) 
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where we have included a constant shift [|| a in the inte- 
gral in Eq. (0) , which does not affect the equality. Eq. (||) 
can be re-written as 



J P{a) W(a, cf>) B(a - a) da, 



where 



W{a,cp) 



(* T \<i>'(<r-d-)) 
(*t|0) 



(6) 



(7) 



The new propagator in Eq. (g) defines a new random 
walk. In each step the walker |</>) is propagated to \<j>') 
by B(a — a): \<p'{a — a)) = B{a — d)\<fi), where a is 
again sampled from P{a). W(a,<p) is a c-number which 
can be accounted for by having every walker carry an 
overall weight factor and updating them according to: 
w<j)i = W(a,<p)w^. Formally the MC representation of 
\^g) in the new random walk is: 



W) 



(8) 



For any choice of the shift a, the new random walk is 
an exact procedure to realize the imaginary time propa- 
gation, in the sense of Eq. (ft). The optimal choice of a 
is determined by minimizing the fluctuation of W(a, (f>) 
with respect to a. To 0(y/r) this yields 



(9) 



With this choice the leading er-dependent term in W is 
reduced to C(r) and, by expanding B(a — a) in |</>') in 
Eq. (Q), we can manipulate W into the following form: 



W(a, 4>) = exp 
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exp[-TE L (</>)], (10) 



where the term El parallels the local energy in real-space 
QMC methods. Both El and the shift a in Eq. (0) are 
independent of any overall phase factor of \<p). 

The weight of the walker in the new random walk is 
determined by El- In the limit of an exact \*S?t), El is a 
real constant, and the weight of each walker remains real. 
The so-called mixed estimate for the energy is phaseless: 



E G 



(11) 



With a general which is not exact, a natural ap- 

proximation is to replace El in Eq.'s (10) and ([ll]) by 
its real part, ReE^,. We have thus obtained a phaseless 
formalism for the random walk, with real and positive 
weights in Eq.'s (|) and (0). 

Despite this, an additional constraint is still required. 
To illustrate the problem we consider the overlap (^t\<P') 
during the random walk. Let us denote the phase of 



{^?T\<i>'{v — <*)) I T\<t>) by AO, which is in general non- 
zero (of order —alma). This means that, the walkers 
will undergo a random walk in the complex plane defined 



by (* 3 



At large (3 they will therefore populate the 



complex plane symmetrically, independent of their initial 
positions. It is useful to contrast the situation with the 
special case of a real v. For any v the shift a diverges as 
a walker approaches the origin in the complex plane, i.e., 
as (^t\4>') ~ * 0- The effect of the divergence is to move 
the walker away from the origin. With a real v, AO = 
and the random walkers move only on the real axis. If 
they are initialized to have positive overlaps with I^t), 
a will ensure that the overlaps remain positive through- 
out the random walk, much like in fixed-node diffusion 
Monte Carlo (DMC) in real space. Thus in this case 
the phaseless formalism reduces to the constrained path 
Monte Carlo method of Ref. || , and it alone is sufficient 
to control the sign problem. For a complex v, however, 
the random walk is "rotationally invariant" in the com- 
plex plane, and the divergence of a is not enough to pre- 
vent the build-up of a finite density at the origin. Near 
the origin the local energy El diverges, which causes di- 
verging fluctuations in the weights of walkers. To address 
this we make an additional approximation. We project 
the random walk to "one-dimension" and multiply the 
weight of each walker in each step by max{0, cos(A#)}. 
Imposing instead Hb{9t\<I>') > gave similar results, but 
with somewhat larger variance. 

We apply the new method to Si atom, molecule, and 
bulk. The Si 4+ ions are represented by a norm-conserving 
LDA Kleinman-Bylander (KB) non-local pseudopoten- 
tial |0. We use periodic boundary conditions, and a 
plane-wave basis with a kinetic energy cut-off E cut — 
12.25 Ry. The error resulting from E cut was estimated 
through LDA calculations and is smaller than the MC 
statistical errors. The pseudopotential can be applied 
in essentially the same way as in plane-wave-based LDA 
calculations Calculations involving v and the local 
part of the pseudopotential are efficiently handled using 
fast Fourier transforms. The separable KB form of the 
non-local pseudopotential makes its application as effi- 
cient as in LDA plane- wave codes. Our \^>t) is a single 
Slater determinant consisting of LDA orbitals. 

In Table [|, we show results for the atom and molecule. 
Additional calculations with a = 22as supercells show 
that finite-size errors at a — 19as were smaller than the 
MC statistical errors. Our calculated Si2 binding energy 
is in excellent agreement with the experimental value jl8) . 

In the bulk calculations, we use fee supercells consist- 
ing of 2, 16, 54 atoms (5209 plane waves). As Fig. |l| 
shows, the new method leads to a large improvement. 
Results for 16 and 54 atoms are shown in Table O. Our 
calculation for 54 atoms took several days on 20 Compaq 
Alpha 667 MHz processors. For the bulk cohesive energy, 
we first included a correction for the independent-particle 
finite-size error from the LDA results. We then corrected 



4 



TABLE I: Total valence energies of Si and Si2, and binding 
energy of Si2- The Si2 ground state is 3 E~ (electronic config- 
uration 5 | 3 J.). Calculations were done at the experimental 
equilibrium bond length of 4. 244a b, in a cubic supercell with 
a = 19a.B (4945 plane waves). Energies are in eV. Error bars 
are in the last digit and are in parentheses. 





Si 


Si 2 


Si 2 E B 


LDA 


-102.648 


-209.175 


3.879 


QMC 


-103.45(2) 


-210.03(7) 


3.12(8) 


Experiment 






3.21(13) 



TABLE II: Cohesive energy of bulk Si. Calculations are done 
for fee supercells with 16 and 54 atoms, at a cxp = 5.43A. 
QMC result at oo is from 54 atoms and includes two finite-size 
corrections: (i) an independent-particle correction of 0.311 
eV from LDA and (ii) an additional Couloumb correction of 
—0.174 eV from Ref. [^, ^2|. A zero-point energy correction 
of —0.061 eV was also added to the calculated results at oo. 
Energies are in eV per atom. Error bars are in the last digit 
and are in parentheses. 



16 



54 



LDA 
QMC. 



Ex periment 



3.836 
3 79(4) 



4.836 
4 51 (3) 



5.086 
4 59(3) 



sibilities exist for further improvement of the method. 
An improved \^t) will give improved results. The free- 
dom to choose the one-particle basis and the form of HS 
transformation, both of which can impact the quality of 
the results, offers significant opportunities. For periodic 
systems it should be possible to generalize the formalism 
to allow k-point sampling. 

In conclusion, we have described a method for ground- 
state QMC calculations that allows the use of any one- 
particle basis. The method is general and applies to any 
Hamiltonian of the form in Eq. (]]]). It provides an ap- 
proximate way to control the phase problem in all AF- 
based QMC methods, while allowing many of their ad- 
vantages to be retained that lead to their applications 
spanning several areas. We have shown that the method 
gave accurate results for systems from an atom to a large 
supercell, using a simple trial wave function. 

We are very grateful to E. J. Walter for help with 
LDA calculations and with several programming issues. 
We thank P. Kent for sending us his DMC data and 
J. Carlson for useful conversations. This work is sup- 
ported by the NSF under Grant No. DMR-973404 1 and 
the Resea rch Corporation, and by ONR Grant dmr 
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for the remaining Couloumb finite-size error |19| ] using 
the results of Kent et al. [^0|. Our result is again in 
excellent agreement with the experimental value (from 
Ref. (l|]). It also compares very well with the result of a 
recent fixed-node DMC calculation j2l]], which also used 
a 54-atom supercell and gave 4.63(2) eV per atom after 
similar finite-size and zero-point energy corrections. 

Without an exact solution to the sign/phase problem, 
reducing the reliance on trial wave functions is clearly 
of key importance to increasing the predictive power of 
QMC. For continuum electronic systems such as our test 
cases above, fixed-node DMC has often been the most ac- 
curate theoretical method |u| . It is encouraging that the 
new method, using simple LDA trial wave functions, gave 
comparable results to DMC. For similar supercells DMC 
often uses trial wave functions with 30-100 additional pa- 
rameters [13]. Obtaining a good enough \*$>t) is instru- 
mental to a successful DMC calculation, and often con- 
stitutes a substantial effort. The quality of \^>t) controls 
the systematic errors from the fixed-node approximation 
and the variance. It also affects errors due to the locality 
approximation |23[| , which has been employed by most 
DMC calculations with non-local pseudopotentials. In 
the new method the latter approximation is eliminated. 
It remains to be seen whether the present method could 
lead to more accurate results than fixed-node DMC for 
continuum systems. This has been possible in some cases 
[B4j with real AF in the Hubbard model. 

We have presented a general framework. Various pos- 
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